%% Orbital Simulator
% By: John M. Murray

%% Script Init
%
tic
clear all;
close all;

%% User Settings
%
steps=2.5e6/200;			% Number of steps to take
step_size=200;		% Seconds
body.size=50;		% Kg
%body.loc=[405696 0 0]*1e3;	% meters
%body.vel=[0 959.6 0];	% meters/second

body.loc=[405696 0 0]*1e3;	% meters
body.vel=[0 959.6 0];	% meters/second

center.mass=5.9736e24;	% Kg
center.size=10;		% visual size

%% Parameters/Constants
%
G=6.67428e-11; % m^(3) * kg^(-1) * s^(-2)

%% Initial Config
%
h1=figure(1);
axis square;
hold on;
axis([-500e6 500e6 -500e6 500e6])
plot3(0,0,0,'.r','MarkerSize',center.size);
body.fig=plot3(body.loc(1),body.loc(2),body.loc(3),'.b');
body.line=plot3(body.loc(1),body.loc(2),body.loc(3));

%% Loop
%
for i=1:steps
	%% Calculate accelerations
	%
	a=G * center.mass / norm(body.loc)^2;
	accel=a * -body.loc / norm(-body.loc);
	
	%% Calculate new Velocity
	body.vel= body.vel + accel * step_size;
	
	%% Move point
	%
	body.loc=body.loc + body.vel * step_size;
	set(body.fig,'xData',body.loc(1),...
				'yData',body.loc(2),...
				'zData',body.loc(3));
			
	%% Update Line
	%
	
	set(body.line,...
		'xData',cat(2,get(body.line,'xData'),body.loc(1)),...
		'yData',cat(2,get(body.line,'yData'),body.loc(2)),...
		'zData',cat(2,get(body.line,'zData'),body.loc(3)));
		
	
	%% Regraph
	%
	if(mod(i,10)==0)
		title(sprintf('%d%% Done at %.2f steps a second',round(i/steps*100),i/toc));
		drawnow;
	end
	axis equal;
	
end
%% Fin
toc

%% Notes
% F=m1*a
% F=G*(m1*m2)/r^2
% a=G*m2/r^2

%%
% Earth
% 5.9736e24 % kg

%%
% Moon
% 405696 % Apogee
% 363104 % Perigee